function varargout = lombpsd(x_time, y_data, x_freq, sig_level);

% LOMBPSD - use Lomb's method to compute normalized PSDs
%
%    [F, PSD] = LOMBPSD(T, X) uses Lomb's method to compute normalized PSDs
%    from the vector X sampled at times T to obtain statistical measures of
%    PSD features (Numerical Recipes, section 13.8).  The normalized PSD is
%    evaluated over the frequency range provided by this code.
%
%    [F, PSD] = LOMBPSD(T, X, F) evaluates the normalized PSD over the range
%    of frequencies specified by the vector F.
%
%    [F, PSD, CONF] = LOMBPSD(T, X, F, P) provides estimates of the normalized
%    PSD confidence levels that exceed the significance levels specified by P.
%
%    LOMBPSD(...) creates a plot of PSD vs F and displays confidence levels
%    if specified.
%
%    Note that this routine only works for individual data sets.  If X and T
%    are matrices, they are reshaped into vectors and one PSD is calculated.

% By:   S.C. Molitor (smolitor@eng.utoledo.edu)
% Date: December 22, 2000
%
%    Note that there are two implementations of this code in LOMBPSD.M.  The
%    first is straight from equations 13.8.4 and 13.8.5.  The second uses
%    trigonometric indentities to reduce the number of calls to SIN() and
%    COS(), which are slow.

% validate arguments

if ((nargin < 2) | (nargin > 4))
   msgbox('Invalid number of arguments', 'LOMBPSD Error', 'warn');
   return
elseif (~isnumeric(x_time) | isempty(x_time))
   msgbox('T must be a numeric vector', 'LOMBPSD Error', 'warn');
   return
elseif (~isnumeric(y_data) | isempty(y_data))
   msgbox('X must be a numeric vector', 'LOMBPSD Error', 'warn');
   return
elseif (prod(size(x_time)) ~= prod(size(y_data)))
   msgbox('T and X must have the same number of elements', 'LOMBPSD Error', 'warn');
   return
elseif (nargin < 3)
   x_freq = [];
   sig_level = [];
elseif (nargin < 4)
   sig_level = [];
elseif ((min(sig_level) <= 0) | (max(sig_level) >= 1))
   msgbox('Significance levels must be between 0 and 1', 'LOMBPSD Error', 'warn');
   return
end

% convert data & time to 1 x N row vectors (mean = 0)

t = reshape(x_time, 1, prod(size(x_time))) - mean(x_time);
x = reshape(y_data, 1, prod(size(y_data))) - mean(y_data);

% calculate frequency range if not provided (0 to Nyquist)
% convert frequency to M x 1 column vector (rad/sec)
% convert data to M x N array (identical rows)
% create M x N array of frequency * time

if (~isnumeric(x_freq) | isempty(x_freq))
   Fs = 1 / mean(diff(t));
   x_freq = 0.5 * Fs * [1 : length(t)] / length(t);
end
w = 2 * pi * reshape(x_freq, prod(size(x_freq)), 1);
h = ones(size(w)) * x;
wt = w * t;

% Method 1: straightforward evaluation of eqns 13.8.4 & 13.8.5
% has two calls to SIN() and COS() on M x N arrays

% calculate M x 1 product of angular frequency and tau
% create M x N array of frequency * (time - tau)
% calculate sine & cosine of this M x N array
% create M x N array of data * sine or cosine
% sum across columns to create M x 1 PSD vector

%wtau = 0.5 * atan2(sum(sin(2 * wt), 2), sum(cos(2 * wt), 2));
%arg = wt - (wtau * ones(size(t)));
%coswt = cos(arg);
%sinwt = sin(arg);
%hcoswt = h .* coswt;
%hsinwt = h .* sinwt;
%psd = ((sum(hcoswt, 2) .^ 2) ./ sum(coswt .^ 2, 2) + (sum(hsinwt, 2) .^ 2) ./ sum(sinwt .^ 2, 2))' / (2 * s2);


% Method 2: evaluate eqns 13.8.4 & 13.8.5 with trig identities
% has single call to SIN() & COS() on M x N arrays

% calculate sin(wt) & cos(wt)
% calculate product of data and sin(wt) & cos(wt)
% calculate products & squares of sin(wt) & cos(wt)

coswt = cos(wt);
sinwt = sin(wt);
hcoswt = h .* coswt;
hsinwt = h .* sinwt;
cos2wt = coswt .^ 2;
sin2wt = sinwt .^ 2;
cossinwt = coswt .* sinwt;

% sum arrays across columns (creates M x 1 vectors)

shcoswt = sum(hcoswt, 2);
shsinwt = sum(hsinwt, 2);
scos2wt = sum(cos2wt, 2);
ssin2wt = sum(sin2wt, 2);
scossinwt = sum(cossinwt, 2);

% calculate PSD from M x 1 sum vectors
% PSD is transposed back to row vector (1 x M)

wtau = 0.5 * atan2(2 * scossinwt, scos2wt - ssin2wt);
coswtau = cos(wtau);
sinwtau = sin(wtau);
shcos = coswtau .* shcoswt + sinwtau .* shsinwt;
shsin = coswtau .* shsinwt - sinwtau .* shcoswt;
s2cos = (coswtau .^ 2) .* scos2wt + 2 * sinwtau .* coswtau .* scossinwt + (sinwtau .^ 2) .* ssin2wt;
s2sin = (coswtau .^ 2) .* ssin2wt - 2 * sinwtau .* coswtau .* scossinwt + (sinwtau .^ 2) .* scos2wt;
psd = (((shcos .^ 2) ./ s2cos) + ((shsin .^ 2) ./ s2sin))' / (2 * var(y_data));

% evaluate confidence levels if specified
% follows exponential distribution w/ M independent measures
% M approximately equals the number of data points if uniform sampling rate
% calculation uses approximation for small values (prevent underflow)

if (isempty(sig_level))
   conf_level = [];
else
   sig_level = reshape(sig_level, 1, prod(size(sig_level)));
   conf_level = ones(size(sig_level));
   small_value = sig_level < 0.01;
   conf_level(find(~small_value)) = -log(1 - (1 - sig_level(find(~small_value))) .^ (1 / length(x)));
   conf_level(find(small_value)) = -log(sig_level(find(small_value)) / length(x));
end

% assign output or create plot

switch (nargout)
case 0
   if (isempty(conf_level))
      plot(x_freq, psd);
   else
      x_fs = [min(x_freq) max(x_freq)];
      y_cl = [conf_level' conf_level'];
      plot(x_freq, psd, '-', x_fs, y_cl, '--');
   end
case 2
   varargout{1} = x_freq;
   varargout{2} = psd;
case 3
   varargout{1} = x_freq;
   varargout{2} = psd;
   varargout{3} = conf_level;
otherwise
   msgbox('Invalid number of outputs', 'LOMBPSD Error', 'warn');
end
return

